Check that all the directories for the .nc files got made
source_dl <- dir(here("data_raw", "CMIP6"))
source_id <- idx$source_id %>% unique() %>% str_to_lower() %>% str_replace_all("-", "_")
stop_if_not(!any(!source_id %in% source_dl))
Check that all the corresponding .csv files exist
csvs <- list.files(here('data'))
stop_if_not(!any(!paste0(source_id, "_data.csv") %in% csvs))
For this analysis, I only want to use models with pr, tas, hfss, and hfls variables in all 5 scenarios (historical, ssp126, ssp245, ssp370, and ssp585)
Perform necessary calculations to compare PET and SPEI among models and between models and observed.
For PET, I’m using the “energy-only” method proposed by Milly and Dune (2016) eq. 8:
\[ PET = 0.8(R_n - G) \]
Except that in their notes, they estimate \(R_n -G\) as hfls + hfss after converting to units of mm/day using the latent heat of vaporazation of water, given by their eq. 2:
\[ L_v(T) = 2.501 - 0.002361T \] in MJ/kg
For the observed data and the CMIP6 data from the same period, I calculate 3-month SPEI using precipitation and PET.
Because SPEI, the variable of interest, is standardized, we looked for models that captured seasonality of precipitation rather than focusing on how well they estimate the exact amounts of precipitation. We calculated mean precipitation for January through December and then calculated a correlation coefficient between the 12 monthly means from each CMIP model and the 12 observed means. We eliminated models with correlation coefficients less than 0.6 for precipitation. Additionally, SPEI takes evapotranspiration into account, specifically through hfls and hfss. Because these variables were not available for observed data, we used temperature but with a less stringent cutoff of correlations greater than 0.4. Additionally we calculated SPEI for each model and counted the number of droughts (SPEI < -1) in each month. However, SPEI was not used to determine which models remained in our ensemble because we did not have strong expectations that accuracy of past SPEI or drought frequency is a good measure of GCM model skill in predicting future SPEI.
cmip_drought <-
df_list %>%
map_dfr(~.x %>%
filter(time >= ymd("1980-01-01") &
time <= ymd("2015-09-30")) %>%
calc_drought_duration, .id = "source_id")
cmip_drought <-
bind_rows(
calc_drought_duration(xa_spei) %>% add_column(source_id = "observed"),
cmip_drought
) %>%
select(source_id, mean_n_mon, sd_n_mon)
#calculate number of droughts per month and mean drought duration
# number droughts per month
drought_df <-
df_list %>%
map_dfr( ~{
.x %>% filter(time >= ymd("1980-01-01") &
time <= ymd("2015-09-30")) %>%
mutate(drought = spei < -1) %>%
filter(!is.na(drought)) %>%
mutate(month = month(time)) %>%
group_by(month) %>%
summarize(n_droughts = sum(drought))
}, .id = "source_id")
xa_drought_seasonality <-
xa_spei %>%
mutate(drought = spei < -1) %>%
filter(!is.na(drought)) %>%
mutate(month = month(date)) %>%
group_by(month) %>%
summarize(n_droughts = sum(drought)) %>%
add_column(source_id = "observed")
drought_df <- bind_rows(drought_df, xa_drought_seasonality)
#create mini-histograms of numbers of droughts
drought_plots <-
drought_df %>%
group_by(source_id) %>%
nest() %>%
mutate(
drought_season = map(data, ~{
ggplot(.x, aes(x = month, y = n_droughts)) +
geom_col(fill = "orange") +
scale_x_continuous("month", breaks = 1:12,
labels = ~month(.x, label = FALSE)) +
labs(y = "# SPEI < -1") +
theme(text = element_text(size = 45))
})) %>%
dplyr::select(-data)
# drought_plots$drought_season[[1]]
| Comparison of observed data to CMIP6 'historical' output | ||||||
|---|---|---|---|---|---|---|
| Data only from 1980 to 2015 to match observed. | ||||||
| Source | Seasonality (monthly means)1 | Droughts (SPEI < -1) | ||||
| pr cor | precipitation | tas cor | temperature | mean ± SD duration (months) | drought seasonality | |
| observed2 | 1.00 | 1.00 | 3.4±2.8 | |||
| cas_esm2_0 | 0.93 | 0.93 | 2.1±1.6 | |||
| fgoals_f3_l | 0.94 | 0.86 | 2±1.4 | |||
| awi_cm_1_1_mr | 0.80 | 0.80 | 2.4±1.3 | |||
| fgoals_g3 | 0.80 | 0.74 | 2.2±1.8 | |||
| taiesm1 | 0.77 | 0.79 | 3±2.5 | |||
| cmcc_esm2 | 0.72 | 0.77 | 2.2±1.5 | |||
| access_esm1_5 | 0.68 | 0.67 | 2.1±1.9 | |||
| cmcc_cm2_sr5 | 0.60 | 0.71 | 2.5±2.2 | |||
| canesm5 | 0.50 | 0.37 | 3.1±2.5 | |||
| bcc_csm2_mr | 0.14 | 0.63 | 2.2±1.3 | |||
| ec_earth3_veg_lr | 0.25 | 0.32 | 2.6±2 | |||
| iitm_esm | 0.25 | −0.04 | 2.2±1.3 | |||
| access_cm2 | 0.05 | 0.12 | 2±1.2 | |||
| cams_csm1_0 | 0.20 | −0.09 | 2±1.2 | |||
|
1
Red numbers highlight correlations (Pearson's r) < 0.6 for precipitation and < 0.4 for mean temperature.
2
Observed data from Xavier et al. (2016)
|
||||||
Below are plots of all data downloaded from each CMIP6 source.